home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
BIOLOGY
/
GSRC208A.ZIP
/
NEWTON.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-07-16
|
4KB
|
109 lines
#include "copyleft.h"
/*
GEPASI - a simulator of metabolic pathways and other dynamical systems
Copyright (C) 1989, 1992 Pedro Mendes
*/
/*************************************/
/* */
/* steady-state solution by */
/* the damped Newton method */
/* */
/* Zortech C/C++ 3.0 r4 */
/* MICROSOFT C 6.00 */
/* Visual C/C++ 1.0 */
/* QuickC/WIN 1.0 */
/* ULTRIX cc */
/* GNU gcc */
/* */
/* (include here compilers that */
/* compiled GEPASI successfully) */
/* */
/*************************************/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "globals.h" /* global symbols */
#include "globvar.h" /* extern global variable */
#include "matrix.h" /* matrix manipulation */
#include "rates.h" /* rates and jacobian calculation */
#include "datab.h" /* database structures */
#define N_OK 0
#define N_LMIN 1
#define N_NOCONV 2
#define N_JSING 3
int newton( void )
{
register int m, i, j, k;
unsigned long df; /* dumping factor (power of two) */
static double xpss[MAX_MET], /* x at iteration m */
ma[MAX_MET][MAX_MET], /* inverse of jacobian matrix */
mb[MAX_MET][MAX_MET], /* copy of jacobian matrix */
h[MAX_MET], h2[MAX_MET], /* increment and damped increment */
re,nre; /* residual error & new res.err. */
for(j=0;j<totmet;j++) /* make x[i] the first guess */
xpss[j] = xss[j] = x[j];
/* initial values of f(x) and re */
calcrates( xpss ); /* f(x) */
for(j=0, re = (double) 0; j<nmetab;j++) /* ||f(x)||2 residual error */
re += rate[j] * rate[j];
re = sqrt( re );
if ( re < options.hrcz ) return N_OK; /* if already in s.s. job is done!*/
for(m=0;m<newtlim;m++) /* iterate "until satisfied" */
{
if (options.debug) printf("\nnewton() - iteration %2d, ", m);
calcjacob( xpss ); /* f'(x) jacobian */
memcpy( mb, jacob, MAX_MET * MAX_MET * sizeof(double) );
lu_inverse( (double(*)[MAX_MET][MAX_MET])mb,
(double(*)[MAX_MET][MAX_MET])ma,
indmet ); /* inv( f'(x) ) */
multmtxv( ma, rate, h,
indmet, indmet, indmet ); /* h = inv(f'(x)) . f(x) */
nre = 2 * re; /* just to start with nre > re */
for( i=0, df=1; (i<32) && (nre>re) ; i++ )
{
if (options.debug) printf(".");
if (i) df *= 2; /* 2 raised to the i power */
for( j=0; j<indmet; j++ )
h2[j] = h[j]/df;
addvctsm( xpss, h2, -1.0, xss, indmet );/* xss = xpss - h2 */
/* update the conserved moieties */
for( j=indmet; j<nmetab; j++)
{
xss[j] = moiety[j];
for( k=0; k<indmet; k++)
xss[j] -= ld[j][k] * xss[k];
}
calcrates( xss ); /* f(x) */
for(j=0, nre=0.0; j<indmet; j++) /* ||f(x)||2 residual error */
nre += rate[j] * rate[j];
nre = sqrt( nre );
}
re = nre; /* keep nre for next step */
if (i==32)
{
if ( re < options.hrcz ) return N_OK; /* this is a zero, after all! */
else /* stuck in a local minima */
{
if (options.debug) printf(" local minima ");
return N_LMIN;
}
}
if (options.debug) printf(" re=% .4le", re );
for(j=0; j<nmetab; j++)
xpss[j] = xss[j];
if ( re < options.hrcz )
return N_OK;
}
if (options.debug) printf(" no solution ");
return N_NOCONV; /* no convergence */
}